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1 .  INTRODUCTION 


1.0  OBJECTIVE 

The  primary  objective  of  this  study  is  to  develop  the  Finite  Element 
Method  (FEM)  for  electromagnetic  techniques  to  solve  the  problem  of 
transient  scattering  directly  in  the  time  domain.  The  techniques  will 
then  be  applied  to  compute  the  time-dependent  currents  induced  on  curved 
thin  wires  due  to  an  arbitrary  incident  plan  wave  pulse.  Special  cases 
of  this  are  straight  thin  wires  and  Gaussian  pulse. 

1.1  Relevance  of  the  Study 

Transient  electromagnetic  response  of  structures  such  as  strategic 
weapon  systems  and  strategic  command,  communication  and  control  systems  to 
a  nuclear  electromagnetic  pulse  are  of  great  concern  from  the  point  of 
view  of  their  vulnerability  and  survival.  Again,  the  importance  of  tran¬ 
sient  response  cannot  be  overstated  in  radar  target  identification,  elec¬ 
tronic  warfare  and  electronic  countermeasures.  For  example,  the  impulse 
response  can  give  a  useful  characterization  of  each  radar  target  since 
such  a  response  contains  all  necessary  radar  information  in  a  compact  and 
understandable  form.  Strategic  systems  should  be  designed  to  survive  the 
nuclear  transients.  Thus,  an  understanding  of  the  response  becomes  manda¬ 
tory  to  impact  on  and  improve  the  designs  of  systems.  Since  the  systems, 
in  general,  are  complicated  structures,  they  can  only  be  solved  numeri¬ 
cally;  hence  efficient  and  economical  numerical  techniques  are  required. 
The  method  developed  in  this  study  is  expected  to  offer  such  a  tool. 

The  approach  is  based  on  a  unique  technique  for  the  computation  of 
currents  and  fields  on  arbitrary  structures  excited  by  arbitrary  sources. 
This  technique  called  "Finite  Element  Method  for  Electromagnetics  (FEMEM)" 
is  predicated  on  a  variational  principle  governing  the  physics  of  the 
problem  and  an  approximation  procedure  to  carry  out  the  variational 
expression  integration. 
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1 .2  A  General  Discussion  of  Transient  Solutions 

There  are  essentially  two  approaches  for  solving  linear  electro¬ 
magnetic  problems.  One  is  an  indirect  approach  in  which  the  physics  of 
the  problem  is  abstracted  either  by  a  differential  or  by  an  integral  equa¬ 
tion  with  frequency  as  the  variable.  The  equations  are  solved  in  the  fre¬ 
quency  domain  and  then  the  time  domain  solution  is  obtained  by  inverse 
Fourier  transform.  In  the  other  approach,  the  governing  equations  are 
formulated  and  solved  in  the  time  domain. 

In  the  time  domain,  the  problem  can  be  formulated  either  in  terms 
of  differential  or  integral  equation.  From  a  numerical  solution  point  of 
view,  the  integral  equation  approach  offers  definite  advantages  over  the 
differential  equation  approach  with  respect  to  solution  stability  and 
imposing  boundary  conditions. 

1 .3  The  Finite  Element  Method  in  the  Transient  Domain 

The  finite  element  method  has  been  successfully  applied  to  a  host 
of  static  or  steady  state  problems,  including  eigenvalue  problems  through¬ 
out  the  many  engineering  disciplines.  The  extension  of  the  method  to 
transient  problems  may  be  credited  to  Wilson  and  Nickel  1^^^  in  their  study 
of  the  heat  conduction  equation.  Most  of  the  early  papers  in  this  area  are 
concerned  with  solutions  to  the  diffusion  equation  in  one  form  or  another. 
Although  the  wave  equation  has  been  considered  generally  by  Oden^^^  there 
appears  to  be  no  specific  solutions  to  this  equation  for  transient  prob¬ 
lems.  In  general,  three  different  approaches  are  used  in  solving  the 
time  domain  problem  in  conjunction  with  the  FEM.  They  are: 

(1)  In  this  method,  the  transient  solution  is  obtained  by  develop¬ 
ing  a  recurrence  relation  with  the  ordinary  finite  element  equations  for 
the  problem  and  then  time-stepping  progressively.  This  technique  will  be 
further  discussed  later. 

(2)  This  method  depends  on  the  idea  of  incorporating  the  time 
dimension  directly  into  the  finite  element  analysis  as  another  one  of  the 
unknown  nodal  degrees  of  freedom  of  the  systems.  In  this  manner,  time  is 
discretized,  as  well  as  the  spatial  variables.  Here  the  time  span  of 


interest  is  divided  into  finite  elements.  Thus,  the  initial  value  prob¬ 
lem  is  converted  to  a  boundary  value  problem.  Solutions  for  all  intervals 
of  time  are  obtained  simultaneously,  with  nodes  on  each  wire  or  surface 
t  «  constant  defining  the  configuration  of  the  system  at  that  time.  The 
increase  in  problem  size  due  to  the  added  time  dimension  is  a  disadvantage, 

(3)  In  this  approach,  the  solution  is  obtained  by  the  mode  super¬ 
position  method.  This  technique  is  also  known  as  the  normal  mode  method 
or  as  modal  analysis.  The  basis  of  this  method  is  that  the  modal  matrix 
of  the  eigenvalue  problem  can  be  used  to  diagonalize  the  problem  and  thus 
decouple  the  multiple  degrees  of  freedom  problem  to  give  several  one- 
degrees  of  freedom  problem. 

One  advantage  of  the  mode  superposition  method  over  the  direct 
integration  methods  is  that  it  reduces  the  number  of  equations  to  be 
solved.  Since  the  lower  normal  modes  play  a  more  significant  role  in  the 
response  .than  the  higher  modes,  only  the  lower  modes  need  to  be  used. 

This  method  has  the  disadvantage  of  requiring  the  eigenvalue  problem  solu¬ 
tion.  Again,  if  the  number  of  degrees  of  freedom  is  large,  the  eigenvalue 
problem  is  difficult.  Superposition  method  is  applicable  only  to  linear 
problems.  Thus,  it  transpires  that  the  modal  superposition  method  is  less 
general  than  the  other  two  methods  mentioned  earlier. 

However,  it  must  be  mentioned  that  the  advantages  inherent  in  the 
finite  element  formulations  can  be  profitably  used  in  all  three  methods. 
This  report  primarily  concerns  itself  with  the  first  method.  The  subse¬ 
quent  sections  discuss  the  problem  formulation,  FEM  methodology,  code 
development,  numerical  solutions,  discusssions,  and  conclusion. 


2.  MATHEMATICAL  FORMULATION 


2.0  FORMULATION  OF  THE  VARIATIONAL  INTEGRAL 

The  application  of  the  FEM  technique  requires  that  we  select  the 
proper  variational  principle  for  the  posed  problem,  express  the  functional 
Involved  In  terms  of  approximate  assumed  current  distribution  functions 
which  satisfy  the  boundary  conditions  and  minimize  this  functional  to 
obtain  a  set  of  governing  equations  which  Is  then  solved  for  the  unknown 
current  distributions  at  the  nodes. 


z 


Figure  2-1.  Geometry  of  the  Problem 

The  wire  Is  Illuminated  by  a  plane  electromagnetic  pulse  E  (t) 
with  arbitrary  polarization  and  angle  of  incidence. 

Here  the  relations  will  be  developed  for  general  validity.  Then, 
these  will  be  specialized  to  the  problem  at  hand.  Let  (r,t)  be  the 
Induced  current  on  the  perfectly  conducting  structure.  The  boundary 
condition  applicable  at  the  surface  of  a  perfectly  conducting  body  Is 
given  by 

nx?^-nx(?*+?^)»0 

A  ■►t  ■►e 

where  n  Is  unit  normal  to  the  surface  and  E  t  E  and  E  are  the  total 
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incident  and  scattered  fields,  respectively.  This  implies  that  the  tan¬ 
gential  electric  field  is  zero. 


The  variational  form  functional  L(J)  governing  the  physics  of  the 
problem  and  containing  the  quantity  of  interest  J  is  given  by 


L(J)  ■  /  y*  iJ(f,t)  •  K  (r.rjt)  •  iJ  (r;t)  dr  dr' 
-2  y  3”  (r,t)  •  (f,t)  dr 


(2-2) 


where^  and^^  denote  Cauchy  principal  value  integrations  over  the  struc¬ 
ture  and  dr  and  dr'  differential  elements.  The  Kernel  K  (r,r' ,t)  is  a 
complicated  <ntegro-differential  operator  and  is  given  by 


K  (r.r 


•t)  *  4—)  A  A*  *  ^  1  S.T 

r2  7  j  dt  .  t 


(2-3) 


where  s,  s  are  the  unit  tangent  vectors  at  r  and  r 

and  H  ,  e  and  n  are  free  space  permeability,  permittivity  and  imped- 

0  0  X)  ^ 

ance,  respectively;  R  =  r  -  r',  the  vector  distance  between  the  observa- 
tion  point  r  and  the  source  point  r';  v  •  denotes  the  divergence  operation 

P 

in  the  source  coordinates;  t  ■  t  -  —  is  the  retarded  time.  It  is  easily^ 

^  c  ,  i 

seen  that  the  variation  of  L(J)  with  respect  to  0  leads  to  the  time  domain 
electric  field  integral  equation 


S  L(J).  M’  (7.t)  -hf  ]ir  ^ ' 


(2-4) 


f  J(r'.t)  dr'  .  0.  T«t-|  . 

0  *<0 

Now  that  the  variational  form  is  set  for  the  problem,  the  remain¬ 
der  of  the  FEM  technique  is  a  procedure  for  rendering  L(J)  stationary  by 
using  an  expression  for  J. 
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2.1  Solution  of  the  Variational  Integral  Equation  by  the  Finite 
Element  Method 

The  FEM  Is  primarily  a  numerical  procedure  for  solving  complex  prob 
lems.  The  method  was  originally  used  in  the  field  of  structural  mechanics 
but  since  its  roots  belong  In  mathematics  as  a  class  of  approximation  pro¬ 
cedure,  it  can  be  applied  to  a  wider  variety  of  problems  in  other  areas. 

In  the  FEM,  the  region  of  interest  is  divided  into  sub-domains  or  finite 
elements,  with  some  functional  representation  of  the  solution  being 
adopted  over  the  elements  so  that  the  parameters  of  the  representation 
become  unknowns  of  the  problems.  Usually  the  element  parameters  are  the 
nodal  values  and  their  derivatives  at  the  nodes.  Although  the  region  of 
the  problem  is  discretized  into  elements,  the  whole  domain  remains  as  a 
continuum  because  of  the  imposed  restriction  on  the  continuity  across  ele¬ 
ment  interfaces.  The  mathematical  procedure  of  solving  (2-2)  by  the  FEM 
is  discussed  in  the  following  sections. 

2.1.1  Segmentation  in  the  Space  and  Time  Coordinates 

Examination  of  Eq  (2-2)  shows  that  the  source  current  at  r  delayed 
by  a  time  )r  -  T'l/c  is  affecting  the  current  at  the  observation  point  r. 
Because  of  this  retardation  effect,  Eq  (2-2)  can  be  solved  as  an  initial- 

I 

valued  problem  by  using  a  time  marching  procedure.  This  phenomenon  can 
best  be  visualized  by  considering  the  space-time  diagram  as  shown  in 
Figure  2-2. 


Figure  2-2.  The  Space-Time  Diagram 


space- time  point;  the 
equations  and  they  separate 
into  the  space  and  time 
time  as 

-  tp  U(r'  -  r^)  U(f  -  tj) 

(2-5) 

(2-6) 


with  and  as  the  spatial  and  temporal  increments  and  represents 
the  current  value  within  the  space  segment  i  and  time  interval  J.  There¬ 
fore  if  one  postulates  that  the  incident  field  and  all  surface  current  on 

S  are  known  or  equal  to  zero  for  all  time  less  than,  say  t  ,  then 

0 

the  retarded  time  effect  allows  us  to  start  the  solution  at  time  t^  and 
to  view  the  integral  equation  as  an  initial-valued  problem  in  a  "marching 
on"  procedure  in  time. 

2.1.2  The  Subdivision  of  the  Spatial  Region  (A  Generalized  Approach) 

The  region  R  is  subdivided  into  discrete  sub-regions  or  elements, 
each  of  the  same  general  form,  as  shown  in  Figure  2-3,  with  the  boundaries 
of  each  element  being  plane  or  curvilinear  faces,  and  with  the  adjacent 
boundaries  of  any  pair  of  elements  being  coincident.  Commonly  used  ele¬ 
ments  for  surfaces  are  triangular  or  polygonal  form.  At  similar  positions 
in  each  element,  a  number  of  points  are  identified  as  nodes.  They  are 
generally  at  the  vertices  of  the  elements,  and  at  positions  such  as  the 
center  of  an  edge,  the  centroid  of  a  face  or  the  centroid  of  the  element 
vol ume . 

Let  us  denote  the  nodal  values  of  the  solution  4  at  the  p^^  node  as 
4p.  Let  the  number  of  elements  into  which  region  R  is  subdivided  be  N^, 


In  the  space-time  diagram  each  dot  represents  a 
solid  lines  are  the  characteristics  of  the  wave 
the  past  and  the  future.  To  divide  the  current 
coordinates,  we  expand  the  current  in  space  and 


*  ^ 


where  U  (r*  -  r^) 


i-1  j-1 
1 


ri.  f 


if  I?'  -?^|<^ 


otherwise 


U  (f  -  t.) 


’  tnt'-tjisji 

0  otherwise 
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ELEMENTS 


Subdivision  of  the  Region  Into  Finite  Elements 


and  the  total  number  of  nodes  In  R  >  D  B  (Boundary)  be  n^  and  n^.  The 
Total  number  of  nodes  In  a  single  element  be  n^.  Then  the  nodal  values  of 
^  can  be  generally  expressed  as  a  column  vector 


♦i 

♦2 


(♦} 


(2-7) 
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2.1.3  The  Element  Shape  Function 

To  solve  Eq  (2-2)  by  the  FEM,  one  needs  to  define  some  shape  func¬ 
tions  or  Interpolation  functions.  These  functions  allow  us  to  express  the 
solution  ^  at  any  position  of  R  In  terms  of  only  the  nodal  values 
Therefore,  we  assume  that  the  solution  ^  can  be  described  In  functional 
forms,  element  by  element,  across  the  region,  I.e.,  can  be  defined  piece- 
wise  over  the  region.  Within  each  element.  It  will  be  supposed  that  ^  can 


be  described  by  a  linear  combination  of  functions  N.  ,  N,  ,  .  .  .  N.  , 

,  .  ,  .  .  .  .  ♦  ,  thus 


Nj®,  and  nodal  values  ^j®. 


♦  "T*  N  *  ♦  *  +  N,®  a  ®  +  N  ®  ♦  ®  +  .  . 
11  2  2  3  3 


N 


or,  in  matrix  notation 

♦  ■  E  (“1'  "2*  ■ 


k  H  ' 
s  »  1,2, 


N  *)  (♦*} 


(N®)  {♦*} 


e 


^  Ns*  C* 
(2-8) 

(2-9) 

(2-10) 


Note  that  the  superscript  e  Is  used  here  to  Identify  a  particular  element. 

The  shape  functions  (N®)  are  restricted  to  being  functions  of  posi¬ 
tions.  Since  the  true  solution  4  Is  prescribed  as  being  continuous  and 
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Is.a.^ 


with  continuous  derivatives  (up  to  some  order)  across  the  region,  the  piece- 
wise  representation  (2-10)  should  have  the  same  properties.  Therefore  the 
shape  functions  are  restricted  by  the  following  conditions: 

1 .  Nj*  ■  1  at  the  node 

2.  *0  outside  element  e,  with  the  node 

as  one  of  its  nodes 

3.  ■  N(r),  a  position  function  within  the  elements. 

In  choosing  the  shape  function,  one  has  to  pay  attention  to  convergence  in 
the  FEM.  Since  it  is  recognized  that  the  FEM  solution  to  a  problem  with  a 
given  size  of  element  is  necessarily  an  approximation  to  the  exact  solu¬ 
tion,  there  must  be  an  assurance  that  successive  finite  element  solutions 
using  smaller  and  smaller  elements  will  converge  smoothly  to  the  exact 
solution  as  the  element  size  tends  to  a  point.  While  comprehensive  condi¬ 
tions  ensuring  convergence  are  not  yet  known  for  all  types  of  linear  problems, 
there  are  certain  criteria  that  must  be  observed  in  order  to  obtain  conver¬ 
gent  solutions: 

(1 )  Completeness 

This  means  that  the  piecewise  representation  (2-10)  within  the 
element  of  the  variable/derivative  in  a  key  integral  must  be  capable  of 
representing  any  continuous  function  as  the  element  size  decreases. 
Mathematically,  the  piecewise  representation  calls  for  a  complete  set  of 
functions  such  as  a  polynomial  function  with  infinite  number  of  terms. 

However,  in  a  FEM  representation,  only  a  finite  number  of  terms  is  taken. 

But  as  pointed  out  by  Melosh  [5]  and  by  Zienkiewicz  [6],  a  monotonic  con¬ 
vergence  can  still  be  obtained  if  the  number  of  terms  used  in  the  repre¬ 
sentation  allow  the  variable/derivative  up  to  and  including  order  t  to 
take  up  any  constant  value  within  the  element,  where  t  being  the  highest- 
order  derivative  of  the  variable  in  the  variational  functional. 

(2)  Compatibility 

This  means  that  the  representation  of  the  variable/derivative  in  a 
key  Integral  of  (2-4)  must  tend  to  the  same  continuity  as  the  exact  solu¬ 
tion,  across  the  inter-element  boundaries,  as  the  size  decreases  to  a  point. 

If  for  a  given  variational  functional,  the  highest-order  derivative 
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Involved  is  of  order  n,  the  derivatives  of  order  up  to  and  including 
(n-1)  are  known  as  the  principal  derivatives  of  that  variable.  Presuming 
that  the  exact  solution  of  the  dependent  variables  are  continuous  with 
continuous  derivatives  up  to  at  least  order  n.  One  weak  requirement  that 
the  compatibility  criterion  is  satisfied  is  to  require  that  the  variable 
and  their  principal  derivatives  are  continuous  in  the  shape  function 
representation.  This  means  that  the  highest-order  derivative  in  a  key 
integral  will  have  a  representation  that  is  at  worst  piecewise  continuous, 
in  which  case  the  representation  will  tend  to  be  continuous  as  the  element 
reduces  to  a  point.  In  general,  completeness  and  compatibility  are  suffi¬ 
cient  conditions  for  convergence  in  variational  finite  element  methods. 
However,  these  conditions  are  very  strong  and  can  be  relaxed  [7].  In 
practice,  the  shape  functions  will  not  be  an  exact  representation  of  the 
true  solution,  but  an  approximate  one,  and  the  solution  obtained  will  be 
similarly  approximate. 

2,1.4  The  Subdivision  of  the  Functional 

Since  (2-2)  represents  essentially  a  quadratic  function,  we  can 
write  it  as 

♦  *  /  •  •  •  •  “d^  (2-11) 


F(u  ,u  ,u 
12  3 


,  Uj)  -  a 


11 


u^-t-a  uu-t'a  uu 
1  12  1  2  13  1  2 


+  a 


21 


u,  + 

1  22 


+  a 


dd 


(2-12) 


and  D  represents  the  domain  of  integration  which  can  be  a  line,  surface  or 
volume,  and  u^,U2,u^  ...  ,u^  represent  the  solution  ^  and  its  various 

derivatives,  ^  ^  ,  ^  .  .  .  .  In  matrix  notation  (2-11)  becomes 

A  j 

♦  {u}'’’[A]  {u}  dD  (2-13) 

D 


where  [A]  is  a  dxd  matrix  and  (u)  a  dxl  column  vector,  or 


*11 

*12  •  •  •  * 

*ld 

[A]  - 

*21 

• 

*22 

*2d 

*d2  •  •  •  • 

*dd 

* 

(2-14) 
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L“^J 

and  superscript  T  denotes  the  transpose  of  a  matrix.  In  general,  the 
matrix  elements  a^j  are  functions  of  the  position. 

If  •*  is  the  contribution  of  an  element  to  the  total  integration  in 
(2-13),  then  this  equation  can  be  written  as 


^  ^  r  T 

e-1  e-1  *'d 

e 

where  represents  the  domain  of  element  e,  let  us  now  consider  a  typical 
term  u^  in  {u}  r  »  0,  1,  2, . By  definition,  u^.  is  a  spatial 

derivative  of  f,  that  is  u_  ■  il-  =  D_^,  where  represents  a.  spatial 

ac'" 

variable  of  concern.  ^ 


From  (2-10)  we  have 
Thus  within  element  e 


u  *  (N®)  in  element  e. 


•  0^*  «  (D^  N®)  {♦®}  -  (U^®)  (2-17) 

A 

where  (U^  )  represents  the  row  vector  for  the  r-derivative  of  the  shape 
function.  So  applying  (2-17)  for  every  element,  we  obtain 

'  {u}  -  (U)  {/)  (2-18) 


where 


V 


O2  N2' 


D.  N 


d  ”2 


(2-19) 


I 


r 
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Substitution  of  (2-18 


which  shows  that  «  is  now  a  function  of  the  n^  nodal  values  ^  ^ 


2.1.5  The  Stationary  Condition 


In  order  to  solve  (2-20)  we  have  to  invoke  the  variational  princi 
The  condition  that  ♦  is  stationary  is  given  by 


2.1.6  The  Element  Matrix  Equation 


To  get  the  element  matrix  equation  we  have  to  combine  (2-23)  and 

a*® 

(2-20).  Considering  the  terra  for  an  element  e  in  (2-20),  we  get 


a* 

Note  that  a  term  will  be  zero  unless  p  is  one  of  the  element  nodes 

identified  by  1 ,  2,  .  .  .  Ic,  .  .  ,  s.  Also  note  that  the  node  identifiers 
1,  2,  3,  .  .  ,  s  are  not  the  same  as  the  system  node  numbers  which  are 
used  to  represent  the  total  number  of  nodes  in  D.  For  example,  if  the 
triangular  element  e  has  its  three  vertices  identified  in  the  system  node 
numbers  as  7,  9,  and  5,  then  we  can  let  its  node  identifiers  (now  s  ■  3) 
as  1  7,  2  9  and  3  5.  Therefore,  the  only  elements  in  the  column 

vector  that  are  non-zero  are  those  that,  in  terms  of  element  node  identi¬ 
fiers,  are 

«.e 

So  (2-22)  reduces  to 


3» 


3» 


a*' 


3  (r) 


3*/3*i® 
J  34/3^2® 


(2-25) 


Letting 

and  using 


a«/3*. 


[B]  -  [Uf  [A]  [U]  . 


g-^  [Q]  {Y}  -  2[Q]  {Y}  , 


we  obtain  from  (2-24) 


*4 


dD. 


where 


[A^]  is  a  s  X  s  matrix  . 


(2-26) 

(2-27) 

(2-28) 


Since  {4®}  is  constant  with  respect  to  the  integration  we  can  write  (2-28) 
as 


9* 

3{4*} 


CA'«]  {♦«} 


(2-29) 


2-11 


where 


(2-30) 


If  we  substitute  (2-19)  into  (2-24)  and  carry  out  the  mathematics,  we  will 
obtain  for  the  ij^^  element  of  [B®]  as 


*  “2",'  (*2.  "j  *  *22  “2  V  * 
+  .  + 


Note  that  in  (2-31)  the  subscripts  on  the  N®  are  in  terms  of  the  node 
identifiers,  not  system  node  numbers. 

Note  that  the  shape  functions  are  explicitly  defined  functions  of 
spatial  variables.  The  integrand  of  a  particular  term,  say 


/  2(D,  N,')  ID,  , 

could  be  evaluated  as  an  explicit  function  of  x,  y  and  z.  If  a^j  are  con¬ 
stant  coefficients,  the  prescribed  integration  over  the  defined  domain  0^ 
of  the  element  would,  in  consequence,  evaluate  the  term  as  a  scalar.  The 
integration,  if  simple,  could  be  carried  out  analytically.  However,  if 
a^j  are  complex  functions  of  x,  y  and  z,  then  the  integration  would  gener¬ 
ally  require  a  numerical  solution.  Therefore,  the  computational  time 


involved  in  a  problem  depends  very  much  on  whether  a^j  are  simple  or  com¬ 
plex  functions. 


2.1.7  The  Boundary  Condition 

It  is  known  in  boundary-value  problems  that  the  solution  is  not 
unique  unless  it  meets  all  the  required  boundary  conditions.  However,  in 
the  variational  finite  element  methods,  if  the  specified  boundary  condi¬ 
tions  are  natural  boundary  conditions  for  the  problem,  then  it  can  be  shown 
that  the  class  of  admissible  functions  is  not  required  to  satisfy  these. 

In  order  to  illustrate  the  treatment  of  the  boundary  condition  in  the 
matrix  equation  (2-29)  let  us  assume  a  Dirichlet  boundary  condition  such 
that 

♦  *  g(x,y,z)  on  B.  (2-32) 


Using  (2-32),  the  n^  nodal  values  (♦p)0  ^or  the  boundary  nodes  on  B  can  be 
calculated  yielding  n^^  equations  of  the  form 


which  implies  that  if  ♦p  satisfies  the  boundary  condition  and 

a*® 

a  constant  value,  then  ■  o  for  an  element  containing  node 


hence  it  is 
p.  Thus  to 


include  the  B.C.  in  the  element  matrix  equation,  the  simplest  procedure  is 

to  replace  the  row  of  the  matrix  [A^*]  in  (2-29)  by  the  row  matrix  of 
(2-33).  In  other  words,  if  p  is  a  boundary-condition  node,  put  zeros  in 
the  p^^  row  of  the  [A^*]  in  (2-29)  except  for  a  1  in  the  diagonal  position 
and  put  in  the  p^^  row  of  the  driving  vector  the  boundary  value  given  by 
(2-32). 
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3.  THIN  WIRE  SCATTERING 

3.0  APPLICATION  OF  THE  FEM  TO  THE  SCAHERING  OF  AN  ARBITRARY 
PULSE  BY  A  THIN  CONDUCTING  WIRE 

The  geometry  of  the  problem  is  given  in  Figure  3-1.  A  perfectly 
conducting  curved  wire  of  length  L  and  radius  a  is  located  in  free  space 
with  one  of  its  ends  at  the  origin  of  the  Cartesian  coordinates. 


$ 


Figure  3-1.  (a)  Geometry  of  Thin-Wire  Scatterer 

(b)  Subdivision  into  Finite  Elements 

The  wire  is  illuminated  by  an  arbitrary  plane  electromagnetic  pulse,  E^(t). 
As  shown  in  Figure  3-2  the  direction  of  propagation  and  the  polarization 
of  the  incident  pulse  is  defined  by  a  triad  (0,  n)  where  0  and  ^  are 
the  ordinary  angles  in  the  spherical  coordinates  made  by  the  propagation 

A 

vector  k,  and  n  is  the  polarization  angle  between  the  electric  vector  and 
the  plane  of  incidence.  Since  the  wire  is  thin  «  1),  we  can  use  the 
thin-wire  approximation  and  assume  that  the  current  flows  only  along  the 
orientation  of  the  wire.  The  surface  integration  in  (2.2)  now  becomes 
a  linear  integration  along  the  wire  whose  arc  length  is  denoted  by  s,  and 
so  the  variational  equation  (2-2)  for  the  thin-wire  is  reduced  to: 
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7  /  • 


K(s ,s '  ,t)  •  Ig(s '  ,t)  ds 'ds 


I5  (s  ,t)  •  E  \s  ,t)  ds 


(3-3) 


where 


i  V  1  i  '^0  ^  3  .  ^qS’R  J_  1  s-R  ^  S  I  *  f 

K(s.s  .t)  •  ^  3t  r2  3s‘  S  ^  r  < 

■  A  ■•/  \  "I#  %i  ^ (3  ^ ) 


and  R  ■  |r(s)-r' (s)  1=  y(r  -  r')  + 

r(s)  is  position  vector  from  the  origin  to  a  point  with  arc  length  s. 

To  convert  the  integral  equation  we  first  divide  the  curved  wire  into  Nj 
uniform  segments  with  number  of  nodes,  and  then  express  the  current  at 
any  point  lying  inside  a  particular  segment  in  terms  of  a  shape  function 
and  the  nodal  values  of  that  segment.  Thus  from  (2-4)  we  have,  after  drop 
ping  the  subscript  s 


i=l  j=l 


Si.  f 


(1 

'f  |S'  -  S,| 

<  41 

U(?' 

-  ^i)  ” 

“  2 

'0 

otherwi se 

U(t' 

if  If  -  tjl 

if 

'  0 

otherwise 

tj)  U(s-  -  s^.)  U(f  -  tj) 

(3-5) 


where  I.|j  represents  the  shape  function  at  the  i^”  segment  and  the 
time  interval.  For  a  given  time  interval,  say  |t'  -  tj^l  , 

I^j(s'  “  s^«  t'  ■  tj)  is  a  function  of  space  only.  Therefore  we  can  express 
the  shape  funct'^on  L,  in  terms  of  spatial  coordinates. 


3.1  Spatial  Shape  Function 

For  the  curved  thin  wire  problem,  v/e  can  represent  the  shape  function 
by  a  polynomial  of  s.  Thus,  for  |t'  -  tj|  i 

I^j(s'  -  Sj.  f  -  tj)  =  Nj  +  N^s  +  NjS^  +  .  .  .  +  \s^~\  (3-6) 

The  coefficients  N|^'s  are  to  be  determined  by  the  continuity  requirement 
across  element  boundaries.  Since  the  variational  equation  (2-2)  involves 
first  spatial  derivative,  it  is  necessary  to  use  a  polynomial  of  at  least 
second  order  in  order  to  meet  the  completeness  and  compatibility  require¬ 
ments  for  convergence.  Thus  we  let 


I^.jCs'  -  s^,  f  -  tp  «  Nj  +  NjS  +  Njs' 


(3-7) 


external 

node 


internal 

node 


Figure  3-3.  An  Element  with  an  Internal  Node. 

In  order  to  determine  Nj,  Nj  and  Nj  uniquely,  we  have  to  pick  an  internal 
node  in  an  element  as  shown  in  Figure  3-3,  and  require  that 


^2  ■  Nj  +  N2S2  N3S2 


(3-8) 


AaN-|-NS.(.NS 
’3  'l  '2  3  3  3  ’ 

or  -  Nj  +  NjS^  +  ,  i  «  1,  2,  3 

where  is  the  nodal  value  of  current  at  the  i^^  node. 
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We  can  write  (3-8)  in  matrix  notation  as 


1  *1 
1  s, 
1  s. 


After  some  matrix  algebra  manipulation  we  get 


(3-9) 


\h\ 

(s,-s)(VS)  (VS)(VS)]  ^2 

(S3-S2)  (sr^s)  /rs 


(3-16) 


where  lsl  *  (sj-  Sj)  (sj-Sj)  (sj-Sj)  •  * 

From  Eq  (3-7)  through  (3-11)  we  obtain  for  an  element  connecting  the  1^ 
and  the  (i  +  1  )^*’  nodes 


(3-11) 


‘-V  *1  *  (3-., 


(^*^+i)(VS'+P 


where  the  subscript  m  denotes  the  internal  node.  Although  the  internal 
node  can  be  placed  at  any  position  within  the  particular  element,  it  is 
usually  located  at  the  midpoint  of  that  element. 

3.2  Time  Derivative  and  Integration  Interpolation 

Since  the  kernel  of  the  variational  integral  equation  (2-2)  contains 
also  first  time  derivative,  it  is  necessary  to  do  temporal  interpolation 
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over  adjacent  time  intervals.  A  second-order  Lagrangian  interpolation  is 
usually  sufficient.  Thus  we  let 


f-tp  .  U  (f-tj) 


(3-13) 


^.j+i 


&  a  ^  i  th 

^4  .*»A 


where  represents  the  nodal  current  at  the  i  node  at  the  y  time 
interval . 

To  avoid  extrapolation  into  the  future,  we  have  to  interpolate  the 
current  at  the  time  step  backwards  to  the  j-1  and  j-2  time  steps  when 

<0.5  such  that 

f (f-t,.  i)(t'-t,.) 

I^.(s'-s.,  t'-t.)  »  ^ ^  . 


Equations  (3-13)  and  (3-14)  can  be  simply  written  as 


(3-14) 


(3-15) 
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•  j+1  ) 

Oj  »  j-2 

■  j 


is  either  given  by  (3-13)  or  (3-14).  From  (3-13)  and  (3-14)  we  can  also 


derive  temporal  derivative  and  integration  as 


3 

3t' 


“  (‘•-'j) 


‘>n  ♦fn  •  <3-'6 


n=n. 


and 


"2 


/I,J  U  (t'-tj)df  . 


(3-17 


At 


n=n. 


where  and  can  be  obtained  easily  from  (3-13)  and  (3-14). 
3.3  Matrix  Equation 


Substitution  of  all  the  pertinent  equations  as  derived  above  into 
(2-2)  yields  the  time  step  (i.e.,  t  =  vAt). 


"s  "s 


i«l  j-1  As^  i  [ 


.L'l 


ds'ds 


s 

/  E^'^'s  Ui-v* 

1.1  4*1  £,  ' 


(3-18 


for  s,  <  s  <  Sj  <  S'  <  Sj,, 

-  'k'  i  r  • 

Note  that  i  and  j  are  used  to  denote  the  and  elements  while  k  is 
used  to  denote  the  retarded  time  interval.  The  actual  time  interval 
is  denoted  by  v.  The  summation  ^  and  ^  denote  the  summation  process 

^i 

over  the  spatial  and  temporal  interpolations  as  given  by  (3-12)  and  (3-15). 

To  cast  (3-18)  into  a  matrix  equation  we  invoke  the  stationary  prop¬ 
erty  of  ^3-18)  by  differentiating  it  with  respect  to  each  nodal  current  at 
the  time  interval  and  setting  the  resulting  equations  equal  to  zero. 

Thus 


(3-20) 


The  final  form  of  (3-19)  can  be  symbolically  written  as 

[Z]  »  (Ej  1  t^t^)"^  • 

where  [23  denotes  the  system  matrix  whose  coefficients  are  functions  of 
space  and  time.  However  its  time  dependence  is  the  same  for  every  time 
interval  (assumed  uniform).  Therefore,  matrix  inversion  is  required  only 
once.  (F)  denotes  a  known  column  vector  containing  information  from  pre¬ 
vious  computation. 

The  boundary  condition  imposed  here  is 

♦iv  *  0  (3-21 ) 

where  i  =  1  and  i  =  n”  . 


3.4  Computation  of  Tangents,  arc  Lengths  and  Distances  for  Curved  Wires 

The  evaluations  of  the  dot  products  s-s'  and  s.R,  the  arc  length  s 
and  the  distance  vector  R  are  carried  out  as  follows: 

For  a  given  smooth  curve  y(u)  defined  by  (3-1)  the  tangent  vector  at  a 
point  p(us)  on  the  curve  is 


A 
S  = 


9y(u)  +  gv(u)  +  g7(u) 
y/gx(u)2  +  gy(u)2  +  gz(u)2 


and  the  unit  normal  vector  n  is  given  by 


(3-22) 


A 

n 


4 

du  *  ^ 


(3-23) 


where  p  *  radius  of  curvature  and  g(u)  are  given  by  (3-2). 

The  distance  along  the  curve  between  two  points  p(u^)  and  p(u2)  is 

s  l/gx^(u)  +  Qy  ^(u)  +  9^  ^(u)  du  (3-24) 

To  find  the  distance  vector  R  »  r-f'  +  an  we  let  and  be  the 

direction  cosines  for  the  radial  vector  from  the  coordinate  origin  to  a 
point  p(u),  then  the  position  vector  r  is  given  by: 
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(3-25) 


I 


Irl  (ax  X  +  ay  y  +  a^  z) 


and 


.Ixiui  ^  ^fy{u) 


.  * 


f,(u) 

!rl 


(3-26) 


+  f/(u)  +  f^^u) 


*  r 


A  A 


thus 


(rax-r'a^^*)  x  +  (rv’^'^y'^  ^  (ra^-r'a^’)  2  +  an  (3.27) 


y  y 


3.5  Numerical  Integration 

By  using  the  FEM,  the  integration  over  the  entire  wire  is  now 
reduced  to  a  sunmation  of  integration  over  the  individual  elements.  The 
integration  in  each  element  is  carried  out  numerically  by  replacing  the 
integration  by  its  Riemann  sum  with  unit-weighting  coefficient.  That  is, 
if  we  divide  the  i^^  element  into  N  subdivisions,  we  have 

N 

f(s)ds  «  23f(Sj)  AS  (3-28) 

j-1 

where  ■  the  domain  of  the  i^*’  element 

As  •  the  size  of  a  subdivision 

Sj  *  the  s  coordinates  of  the  center  of  the 
subdivision  of  the  i^*’  element. 

3.6  Matrix  Inversion 

Since  the  problem  is  solved  as  an  initial-value  problem,  it  is  not 
necessary  to  invert  the  matrix  at  each  time  step  of  solution.  Matrix 
inversion  is  done  only  once  at  the  first  time  step  and  the  inverted  matrix 
is  stored  to  be  used  for  the  following-on  time  steps.  Thus  the  solution 
after  the  first  time  step  can  be  written  as 

(♦,,)  ■  I  ft,)  ♦  "'>] 
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4.  COMPUTER  PROGRAM 


4.0  A  BRIEF  DESCRIPTION  OF  THE  COMPUTER  PROGRAM 

The  computer  program  developed  for  this  study  is  called  "TWFEM". 

It  is  written  in  the  Fortran  IV  language.  The  program  consists  of  a  main 
program  and  ten  subroutines.  The  input  to  the  programs  are:  the  maximum 
value  of  the  parametric  variable  for  the  wire,  the  radius  of  the  wire,  the 
number  of  elements  into  which  the  wire  is  divided,  the  size  of  the  time 
step,  the  final  time  for  the  run,  polarization  and  incident  angles,  the 
kind  of  incident  pulse,  pulse  parameters,  and  a  few  control  option 
parameters  for  running  the  program.  The  output  of  the  program  is  the 
current  distribution,  the  indident  field  strength  and  the  segment  excita¬ 
tion  on  each  node  at  each  time  step.  Most  of  these  outputs  are  stored  in 
a  magnetic  tape  and  can  be  saved  for  future  use.  Because  of  this,  the 
program  can  use  the  results  of  the  final  time  step  in  a  previous  run  as 
the  initial  values  for  the  new  run.  This  capability  is  designed  to  save 
computational  time  by  eliminating  duplicated  computation.  The  numerical 
integration  is  performed  by  a  simple  trapezoidal  quadrature,  and  the 
matrix  inversion  is  done  only  once  using  the  gaussian  elimination 
algorithm.  To  save  computational  time,  many  parameters  are  stored  in 
common  blocks. 

4.1  Flow  Chart 

The  structure  of  the  computer  program  is  given  in  a  flow  chart  as 
shown  in  Figure  4-1.  A  sample  print-out  is  given  in  Figure  4-2. 
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Figure  4-2.  Sample  of  print  out. 
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5.  RESULTS  AND  DISCUSSIONS 


5.0  TYPES  OF  INCIDENT  PULSES 

In  this  study  we  consider  four  different  types  of  incident  pulses. 
They  are  defined  as  follows: 


(a)  Gaussian  pulse 
E’(t)  » 


(5-1) 


where  p  *  spread  parameter  (sec"^ ) 

tmax  *  which  the  pulse  reaches  maximum  value 

(sec) 

(b)  Double  exponential  pulse 

E’(t)  =  e"“^  -  e"^^  t  >  0 

where  a  *  4.0  x  10®  (sec"^) 

3  =  4.76  X  10®  (sec'^ ) 

(c)  Rectangular  pulse 

E^(t)  »  1  0  <  t  <  tp 

=  0  otherwise 

where  tp  *  pulse  width  (sec) 

If  tp  is  large,  it  becomes  a  unit  step  pulse. 

(d)  Ramp,  pulse 

E^(t)  -  t  0  <  t  <  t^ 

-  1  t^  <  t 

where  t^.  *  rise  time,  (sec) 


(5-2) 


(5-3) 


(5-4) 
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?.l  TRANSIENT  CURRENTS 

5.1.1  Semi-Circular  Wire  (x  =  0,  y  =  0.159  cosu,  z  =  0.159  sinu) 

Figures  5-1  to  5-6  present  the  induced  currents  as  a  function  of 

time  at  the  middle  of  a  thin  semi-circular  wire  illuminated  by  different 

types  of  pulses  at  various  angles  of  incidence.  The  parameter  of  the 

pulses  are  p  =  9  x  108  sec"^ ,  t^ax  “  ^  nsec,  tp  =  10  nsec  and  t^  =  10  nsec. 

The  length  of  the  wire  is  L  =  0.5  M  and  the  radius  is  given  by  the  shape 

factor  n  =  Zln  (— )  =10.  The  wire  is  uniformly  divided  into  5  elements 
d 

with  eleven  nodes  (six  external  and  five  internal).  The  time  step  is 
it  =  0.167  nsec  which  is  approximately  equal  to  is/c.  The  current  is 
defined  positive  when  it  flow  from  s  *  o  to  s  =  L.  Examinations  of  the 
plots  reveal  many  interesting  points  of  physics  concerning  the  transient 
response  of  thin  wire  structures.  For  example: 

(a)  Even  for  curved  wires,  the  current  displays  damped  oscillations 
at  a  dominant  frequency  which  is  close  to  the  lowest  frequency 
of  a  straight  wire  of  the  same  length. 

(b)  In  general,  the  temporal  development  of  the  current  along  the 
wire  is  a  very  complicated  thing.  It  depends  on  many  factors, 
such  as  the  wire  length,  the  type  of  incident  pulse,  the 
angle  of  incidence  as  well  as  the  shape  of  the  wire.  The 
build  up  of  the  current  can  be  see  as  follows.  First,  the 
current  starts  to  build  up  from  the  end  of  the  wire  where  the 
incident  field  pulse  hits  initially.  As  time  goes  on,  the 
other  part  of  the  wire  is  also  illuminated  and  the  current 
pulse  begins  to  travel  with  the  velocity  of  propagation  to¬ 
wards  the  other  end.  When  the  current  pulse  reaches  the  other 
end  the  current  pulse  reverses  its  direction  of  propagation  as 
the  current  cannot  flow  forward  any  farther.  This  phenomenon 
goes  on  and  on  until  the  current  completely  decays  due  to 
radiation  loss.  At  any  instant  of  time  the  current  is  a 
superposition  of  response  due  to  direct  excitation,  reflections 
along  wire  and  scattering  from  other  parts  of  wire. 

(c)  As  shown  in  Figures  5-2  and  5-3,  it  is  observed  that  any  sudden 
change  in  the  incident  pulse  would  induce  some  new  charges  and 
in  turn  alter  the  current  distribution  on  the  wire.  In  the 
case  of  the  rectangular  pulse  whose  pulse  width  is  10  nsec,  this 
sudden  change  occurs  at  the  end  of  the  pulse  (ct/L  »  6). 
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(d)  The  reason  ^that  the  second  current  peak  is  larger  than  the  first 
peak  in  the  gaussian.  pulse  is  attributed  to  its  relatively  large 
spread  parameter  which  produces  a  very  rapid  rise  and  fall -off 
behavior. 

(e)  For  non-symmetrical  excitation  as  shown  in  Figures  5-4  to  5-6, 
the  reflections  from  both  ends  are  not  equal  and  also  arrive 
at  an  observation  point  at  different  times.  Therefore,  the 
current  development  is  no  longer  dominated  by  the  fundamental 
resonant  frequency  but  contains  higher  harmonics  as  well. 

5.1.2  Parabolic  Wire  (x  =  o,  y  =  0.5u,  z  *  0.25u^) 

Figures  5-7  to  5-12  show  the  transient  currents  at  different  positions 
along  a  parabolic  wire  of  length  L  =  0.574u  and  shape  factor  -  10.01 
under  illumination  by  a  rectangular  pulse  (tp  =  5nsec)  and  a  double 
exponential  pulse.  The  starting  time  is  taken  to  be  zero  when  the 
pulse  first  hits  the  higher  end  of  the  wire  (s  =  L).  As  expected  the 
response  of  the  rectangular  pulse  has  more  ripples  and  roughness  due  to 
its  discontinuities. 

5.1.3  Helical  Wire  (x  =  0.25  cosu,  y  =  0.25  sinu,  z  =  0.25u) 

Figures  5-13  to  5-16  present  the  transient  currents  at  different 
locations  on  a  helical  wire  of  length  L  =  l.lllM  under  a  gaussian  pulse 
(p  =  9  x  108  sec,  tniax  =  1  nsec).  Again  it  is  seen  that  currents  display 
a  more  symmetrical  oscillation  pattern  for  points  at  region  near  the 
middle  of  the  wire,  where  reflections  from  both  ends  are  more  equal, 
than  those  at  positions  close  to  the  ends  of  the  wire. 

5.1.4  Straignt  Wires  (x  =  0,  y  =  0,  z  =  u) 

Figures  5-17  to  5-19  present  transient  currents  on  a  thin  straight 
wire  illuminated  by  a  gaussian  pulse  at  three  different  angles  of  incidence. 
The  wire  is  along  the  z  axis  and  the  incident  wave  is  in  the  y-z  plane 
($  a  90®  and  n  *  0).  It  is  seen  that  the  current  exhibits  strong 
oscillation  at  the  lowest-characteristic  frequency  of  the  wire.  The 
comparison  between  the  FEM  results  and  those  obtained  by  using  the 
method  of  moments  code^®^  are  also  included  In  Figures  5-18  and  5-19. 
Extremely  good  agreement  Is  obtained  between  these  two  different  approaches. 
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5.2  CURRENT  DISTRIBUTION  ON  THE  ENTIRE  WIRE 


Figures  5-20  through  5-85  present  the  snap  shot  type  current 
distributions  on  the  straight  wire,  semi-circular,  parabolic  and  helical 
wires  for  various  types  of  incident  pulses  and  angles  of  incidence  as 
discussed  earlier  in  the  previous  section.  These  plots  clearly  indicate 
the  build-up  mechanism  for  the  transient  currents  on  wire  structures 
at  each  instant  of  time,  say  t  =  T,  the  current  at  a  point  on  this  wire 
is  a  complicated  combination  of  three  different  excitation  mechanisms 
evaluated  in  the  appropriate  retarded  time  frame.  One  source  of 
excitation  is  of  course  the  direct  incident  field,  the  second  is  due 
to  reflections  at  ends  of  the  wire  and  the  third  comes  from  scattering 
from  other  part  of  the  wire.  The  field  near  the  time  of  arrival  of  the 
incident  wavefront  is  determined  by  its  high-frequency  content.  However, 
at  later  times,  long  after  the  wavefront  has  traversed  the  scatterer,  the 
induced  currents  set  up  oscillations  at  the  natural  frequencies  of  the 
wire.  Since  the  wire  is  in  free  space,  leakage  of  energy  to  infinity 
leads  to  damped  oscillations,  and  the  most  weakly  damped  (probably  the 
lowest  mode)  dominates  the  late  time  response.  This  implies  that  the 
long-time  response  of  the  wire  is  determined  primarily  by  its  overall 
size  rather  than  its  detailed  shape.  The  opposite  is  true  for  the 
early-time  response. 

5.3  CONVERGENCE  AND  COMPUTATIONAL  TIME 

Figure  5-86  shows  the  convergence  test  for  the  current  at  the  middle 
of  a  semi-circular  thin  wire  of  length  L  ®  0.5  m  and  o  =  10.  The  number 
of  elements  used  in  the  test  are  4,  5  and  6  which  correspond  to  9,  11  and 
13  nodes.  It  is  seen  that  the  numerical  results  show  a  better  convergence 
in  early  time  than  in  late  time.  This  is  understandable  because  at  late 
time  the  incident  pulse  has  more  or  less  died  down  and  the  driving  source 
is  merely  due  to  the  smaller  residual  current  and  charge  flowing  along 
the  wire  which  are  more  sensitive  to  element  size.  -Also  in  an  iterated 
numerical  solution,  errors' propagate  and  accumulate  as  time  goes  on. 

The  step  sizes  used  in  this  convergence  test  are  chosen  such  that 
*  1,  where  At  is  the  time  step  and  as  is  the  element  size.  The 
computational  time  for  each  run  depends  on  the  number  of  elements  and 


the  number  of  time  steps  used.  A  typical  run  in  this  study  uses  five 
elements  (11  nodes)  and  150  time  steps  and  it  takes  about  100  sec  on  the 
CDC  6600  machine.  However,  it  must  be  noted  that  the  computer  program 
has  not  been  optimized  to  take  into  account  various  factors  such  as 
structure  symmetry  for  the  broadside  illumination  and  possible  analytical 
integration.  Once  this  is  done,  the  computational  time  can  be  reduced 
further. 


CURRENT  (MR) 


rj  rn  co 

CD  GT  tn 


Figure  TRANSIENT  CURRENT-RECTANGULAR  PULSE 


CURRENT  (Nfl) 


figar#  ^-3.  TRANSIENT  CURRENT-UNIT  STEP  PULSE 


»<> 


n 


CURRENT  (m) 


3Flgo»»65-8.  TRANSIENT  CURRENT-RECTRNGULflR  PULSE 


5-13 


RENT-RECTnNGULRR  PULSE 


i-14 


CURRF:NT  (MR) 


5-16 


CURRENT  (MR) 


L  -  0.574 
S  -  0.517  M 

n  -  10.01 

AT-  0.192  NSEC 
e  -  0.000  DEG 
(p  “  0.000  Dl-G 
r)  -  90.00  DEG 


TRRNSIENT  CURRENT-DOUBLE  EXPONENTIAL  PULSE 


CURRENT  (MR) 


CURRENT  (MR) 


Ftiui?e$-V^  TRANSIENT  CURRENT-GflUSSIRN  PULSE 


5-19 


CURRENT  (MR) 


Figure  5-16.  Transient  Current  -  Gaussian  Pulse 
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Figure  5-18.  Transient  Current  -  Gaussian  Pulse 
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Figure  5-31.  Current  Distribution  -  Gaussian  Pulse 
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Figure  S-69.  CURRENT  DISTRIBUTION-RECTANGULAR  PULSE 


5-74 


CURRENT  (MR) 


3.0 

1 

0.2 

1 

0.4 

1 - 1 - 1 

0.6  0.8  1.0 

1 

- 

S  (M) 

L  -  1.111 

M 

T  -  0.186 

NSEC 

n  -  10.01 

AT-  0.371 

NSEC 

e  -  90.00 

DEG 

4)  -  45.00 

DEG 

T)  -  90.00 

DEG 

Fifluf#  5-72.  CtJRKEWT  OfST^fettTi-Of^-GWSSmN  fULSC 


CURRENT  (m) 


CURRENT  (Mfl) 


I 

Figur#  5-74.  CURRENT  OISTRIBUTION-GftUSSIRN  PULSE  I 


5-79 


/SRSPWiT 


CURRENT  (Mfl) 


L  -  l.lll  M 
T  -  4.637  NSEC 

n  -  10.01 

AT-  0.371  NSEC 
e  -  90.00  OEG 
®  -  45.00  DEG 
r,  -  90.00  DEG 


Fi8urtV76.  CURRENT  DISTRIBUTION-GRUSSIRN  PULSE 


CURRENT  (MR) 


5-82 


CURRENT  (MR) 


CURRENT  (MR) 


Figure  5-79.  CURRENT  DISTRIBUTION-GRUSSIRN  PULSE 


5-84 


r  t 


CURRENT  (MR) 


CURRENT  (MR) 


5-86 


CURRENT  (NR) 


Plgurt  S-32;.  CURRENT  DISTRIBUTION-GRUSSIRN  PULSE 


CURRENT  (MR) 


CURRENT  (MR) 


Figure  «-84.  CURRENT  DISTRIBUTION-GRUSSIflN  PULSE 


S-89 


CURRENT  (MFD 


5.0 

0.2 

0.4 

0,6 

0.8 

1.0 

1 

5  (M) 

L  - 

1.111 

M 

T  - 

15.77 

NSEC 

n  - 

10.01 

AT- 

0.371 

NSEC 

e  • 

90.00 

DEG 

9  - 

45.00 

OEG 

•n  • 

90.00 

DEG 

Figurt'5-45.  CURRENT  DISTRIBUTION-GRUSSIflN  PULSE 


5-90 


CURRENT  (MR) 


6:  CONCLUSION 

In  this  study  the  finite  element  method  for  electromagnetic  tech¬ 
niques  has  been  developed  to  solve  transient  scattering  problems  directly 
In  the  time  domain.  The  problem  Is  formulated  In  terms  of  a  variational 
time-dependent  Integro-differentlal  equation  which  Is  to  be  solved  by  a 
finite  difference  scheme  In  time  and  a  finite  element  technique  In  space. 
Based  on  this  approach  a  computer  program  Is  written  to  calculate  the 
transient  current  on  curved  thin-wire  scatterers  when  excited  by  an  arbitrary 
plane  wave  pulse.  Numerical  results  show  good  accuracy  and  convergence  for 
the  FEM  approach.  Thus,  It  transpires  that  the  FEM  can  be  a  good  numeri¬ 
cal  tool  In  solving  transient  electromagnetic  problems.  As  a  numerical 
method  for  solving  electromagnetic  scattering  problems,  the  FEM  offers  the 
following  advantages: 

(1)  Since  the  formulation  Is  based  on  the  variational  princi¬ 
ple,  the  solution  Is  more  stable  and  the  error  Is  minimized. 

(2)  Although  the  region  Is  divided  Into  finite  elements,  the 
whole  domain  remains  as  a  continuum  because  of  the  Imposed 
restriction  on  the  continuity  across  element  Interfaces. 

This  Is  contrary  to  the  point-matching  solution  used  In 
the  method  of  moments  where  the  true  solution  Is  valid 
only  at  the  matching  points  In  the  whole  domain. 

(3)  FEM  approach  Is  particularly  useful  in  handling  complex 
geometries. 

In  spite  of  Its  advantages,  the  FEM-based  time-domain  code  developed 
here  has  two  shortcomings  from  a  numerical  solution  point  of  view.  The 
shortcomings  are: 

(1)  The  mathematical  and  bookkeeping  aspects  of  the  FEM  are 
Involved,  and 

(2)  The  computational  time  seems  to  be  longer  as  the  code  Is 
not  optimized. 

It  Is  therefore  hoped  that  further  research  work  In  this 
area  would  alleviate  these  difficulties. 
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